学习如何将事件以TTree的格式存储到ROOT文件
束流入射到靶上产生核反应,靶上生成的中子和$\gamma$入射到长条型塑料闪烁体探测器的不同位置,并在不同深度与探测介质产生相互作用。
束流的时间信息$T_{Beam}$由Timing Scintillator(薄塑料闪烁体+PMT)给出。
入射粒子形成的光信号向探测器两端传输,经PMT放大形成输出信号。
输出信号分成两路,一路经CFD甄别后形成时间信号 tu 和 td ,送入TDC的stop端;另一路延迟送入QDC,得到能量信号 Qu,Qd 。
获取系统的触发为 $T_u\times T_d\times T_{Beam}$ , 触发信号中$T_{Beam}$的时间为主时序。
触发信号经适当调节送到TDC的Start端和QDC的Gate端。
将入射粒子的原始信息,以及探测器信息逐事件存入root文件的TTree结构中.
入射粒子原始信息:粒子种类、能量,粒子在探测器上的入射位置、入射深度等
探测器信息: 两端时间信息 tu 、 td ,能量沉积信息 Qu , Qd
打开root文件,进行基本的分析。
闪烁体探测器长度:2L
光在闪烁体的传播速度: $v_{sc}$,光衰减常数:$\lambda$
$d=\sqrt{D^2+x^2}$
中子飞行时间:$TOF_n(ns)=\frac{72 \cdot d(m)}{\sqrt{E_n(MeV)}}$, gamma飞行时间:$TOF_{\gamma}(ns)=3\cdot d(m)$
探测器两端时间:$t_u=TOF+(L-x)/v_{sc}$, $t_d=TOF+(L+x)/v_{sc}$
探测器两端能量(x 处能量沉积Q):$Q_u=Q_0 exp[-(L-x)/\lambda]$, $Q_d=Q_0 exp[-(L+x)/\lambda]$

1. 声明所需的变量;
2. 建立TTree的存储;
3. 建立具体的分支;
4. 单个时间测存储;
5. 循环时间的存储;
6. 修改为root脚本运行。
//声明需要注意的问题:
//1、“同样可以使用const double aa=500.;进行声明”,其中加分号终端不会输出;
//2、100后不加点为整型输出,加小数点为浮点型。
//3、可以在后面一直加声明,不需要再阐述数据类型,但是不建议太多,容易乱
const Double_t Dis=500.,Len=100.,del_Dis=5.; //单位为cm,
const Double_t Lambda=380.; //cm, 光的衰减常数
const Double_t V_sc=7.5; //ns/cm, 光在闪烁体中的传播速度
const Double_t Rg=0.3; //出射gamma射线和中子的比值 //假设为均匀分布
const Double_t En0=100.; //MeV; 中子初始能量
const Double_t Eg0=1; //MeV; γ射线的能量
const Double_t EnRes=50.;//MeV
const Double_t TiRes=1.; //ns
const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin.
Double_t Tu_off=5.5; //time offset,//束流飞行时间+PMT的渡越时间+电缆上的传输时间
Double_t Td_off=20.4; //time offset
//声明所需要的参数,即需要输出的参数
//* 入射粒子原始信息:粒子种类、能量,粒子在探测器上的入射位置、入射深度等
//* 探测器信息: 两端时间信息 tu 、 td ,能量沉积信息 Qu , Qd
Double_t X,Energy; //应该是均匀分布
Int_t PID; //pid=1 is neutron、pid=0 is gamma
Double_t T_u,T_d,Q_u,Q_d;
Double_t TOF,cTOF; //TOF:粒子实际飞行时间,cTOF:计算得到的TOF
//声明中间参量,不需要输出
Double_t Deepth; //粒子的飞行深度
Double_t Dis_Deepth;
Double_t Range;
//新建tree的文件,然后声明tree
TFile *Output_File=new TFile("tree.root","recreate");
TTree *Output_Tree=new TTree("tree","tree structure");
Output_Tree->Branch("PID",&PID,"PID/I");
Output_Tree->Branch("X",&X,"X/D");
Output_Tree->Branch("Energy",&Energy,"Energy/D");
Output_Tree->Branch("T_u",&T_u,"T_u/D");
Output_Tree->Branch("T_d",&T_d,"T_d/D");
Output_Tree->Branch("Q_u",&Q_u,"Q_u/D");
Output_Tree->Branch("Q_d",&Q_d,"Q_d/D");
Output_Tree->Branch("TOF",&TOF,"TOF/D");
Output_Tree->Branch("cTOF",&cTOF,"cTOF/D");
TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
TRandom3 *gr=new TRandom3(0);
//开启循环
for(Int_t i_event=0;i_event<100000;i_event++){
X=gr->Uniform(-Len,Len);
Deepth=gr->Uniform(-del_Dis/2,del_Dis/2);
Dis_Deepth=Dis+Deepth;
Range=TMath::Sqrt(X*X+Dis_Deepth*Dis_Deepth);
//判断为哪种粒子
if(gr->Uniform()<Rg){
PID=0;
Energy=Eg0;
TOF=10/3*(Range*0.01);
}
else{
PID=1;
Energy=gr->Gaus(En0,EnRes/2.355);
TOF=72.29824/TMath::Sqrt(Energy)*(Range*0.01);
}
//该给出时间信号了
T_u=TOF+(Len-X)/V_sc+gr->Gaus(0,TiRes/2.35)+Tu_off;
T_d=TOF+(Len+X)/V_sc+gr->Gaus(0,TiRes/2.35)+Td_off;
cTOF=(T_u+T_d)/2.;//simplified calculation.
hctof->Fill(cTOF);
//该存储能量信号了
Double_t q0=Energy*gr->Uniform();
Q_u=q0*TMath::Exp(-(Len-X)/Lambda);
Q_u=gr->Gaus(Q_u,Q_u*QRes/2.35);//QRes, relative energy resolution
Q_d=q0*TMath::Exp(-(Len+X)/Lambda);
Q_d=gr->Gaus(Q_d,Q_d*QRes/2.35);
Output_Tree->Fill();
}
hctof->Write();//写入预定义的histogram到文件
Output_Tree->Write();//写入TTree到文件
Output_File->Close();//关闭文件
%jsroot on
TFile *Input_File=new TFile("tree.root");
Input_File->ls()
tree->Print()
tree->Show(0);
tree->Show(3)
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("T_d-T_u>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
tree->Draw("TOF>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
TH1D *hh=(TH1D*) Input_File->Get("hctof"); //在代码中得到文件内数据指针的方法,在ROOT环境下直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("cTOF:T_d-T_u>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw()
tree->Draw("cTOF:T_d-T_u>>(2000,-20,50,50,40,45)","PID==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯?
tree->Draw("T_d-T_u:X","","colz");//观察两个参量的关联
c1->Draw();
tree->Draw("log(Q_u/Q_d):X","","colz");
c1->Draw()
!jupyter nbconvert 1.1_Creat_TTree.ipynb --to html